# TODO: Add comment
# 
# Author: yaping
# Time: 2013-6-20 
###############################################################################
#give an input file, point out which column it contains a p value, return another file whose last column is the FDR corrected P value, while the other part is the same from the input file.
#by default, the column contain the p value is the last column

source("/home/uec-00/yapingli/code/mytools/R/heatmapRUtils.R")
source("/home/uec-00/yapingli/code/mytools/R/generateHeatmapForNOMeSeq.R")
library(timeSeries)
#default value
#jet.rev.colors <-c("#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F")
jet.rev.colors <-c("#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000","#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000","#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000","#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000","#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F")
jet.colors <-rev(jet.rev.colors)

files<-NULL
heatMapCols<-NULL
ylabForAvePlot<-NULL
rowSideFiles<-NULL
RowSideColorsName<-NULL
colToUses<-NULL
keys<-NULL
capLimit<-NULL
logScale<-NULL
autoScale<-NULL
clusterNum<-NULL

prefix="heatmap"
scale=1000
bin_size_align=1
move_step=20
heatmap_clustering_scale = 1000
heatmap_clustering_bin_size = 20
heatmap_clustering_bin_size_align = 1
heatmap_row_order<-NULL

fileNumToPrintSideBar=1
outputSubClusterCordinate=FALSE
regionToClusterLow=-500
regionToClusterHigh=500
pdfOut=TRUE
twoStepClustering=FALSE
preClustering=FALSE
orderByMaxOccDisToTss=FALSE
multiSampleClustering<-NULL

y_max<-NULL
y_min<-NULL
breaks<-NULL 
wd="./"
prefix<-NULL
for (e in commandArgs(TRUE)) {
	ta = strsplit(e,"=",fixed=TRUE)
	if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "wd"){
			wd<-ta[[1]][2]
		}
		if(ta[[1]][1] == "inputFn"){
			files<-c(files,ta[[1]][2])
		}
		if(ta[[1]][1] == "prefix"){
			prefix<-c(prefix,ta[[1]][2])
		}
		if(ta[[1]][1] == "clusterNum"){
			clusterNum<-c(clusterNum,as.numeric(ta[[1]][2]))
		}
		if(ta[[1]][1] == "heatMapCols"){
			heatMapCols<-c(heatMapCols,ta[[1]][2])
		}
		if(ta[[1]][1] == "rowSideFiles"){
			rowSideFiles<-c(rowSideFiles,ta[[1]][2])
		}
		if(ta[[1]][1] == "colToUses"){
			colToUses<-c(colToUses,ta[[1]][2])
		}
		if(ta[[1]][1] == "RowSideColorsName"){
			RowSideColorsName<-c(RowSideColorsName,ta[[1]][2])
		}
		if(ta[[1]][1] == "keys"){
			keys<-c(keys,as.logical(ta[[1]][2]))
		}
		if(ta[[1]][1] == "ylabForAvePlot"){
			ylabForAvePlot<-c(ylabForAvePlot,ta[[1]][2])
		}
		if(ta[[1]][1] == "scale"){
			scale<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "bin_size_align"){
			bin_size_align<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "move_step"){
			move_step<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "logScale"){
			logScale<-c(logScale, as.logical(ta[[1]][2]))
		}
		if(ta[[1]][1] == "capLimit"){
			capLimit<-c(capLimit, as.logical(ta[[1]][2]))
		}
		if(ta[[1]][1] == "fileNumToPrintSideBar"){
			fileNumToPrintSideBar<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "outputSubClusterCordinate"){
			outputSubClusterCordinate<-as.logical(ta[[1]][2])
		}
		if(ta[[1]][1] == "autoScale"){
			autoScale<-c(autoScale,as.logical(ta[[1]][2]))
		}
		if(ta[[1]][1] == "y_min"){
			y_min<-c(y_min, as.numeric(ta[[1]][2]))
		}
		if(ta[[1]][1] == "y_max"){
			y_max<-c(y_max,as.numeric(ta[[1]][2]))
		}
		if(ta[[1]][1] == "regionToClusterLow"){
			regionToClusterLow<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "regionToClusterHigh"){
			regionToClusterHigh<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "preClustering"){
			preClustering<-as.logical(ta[[1]][2])
		}
		if(ta[[1]][1] == "twoStepClustering"){
			twoStepClustering<-as.logical(ta[[1]][2])
		}
		
		if(ta[[1]][1] == "multiSampleClustering"){
			multiSampleClustering<-c(multiSampleClustering,as.numeric(ta[[1]][2]))
		}
		if(ta[[1]][1] == "pdfOut"){
			pdfOut<-as.logical(ta[[1]][2])
		}
		if(ta[[1]][1] == "breaks"){
			breaks<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "heatmap_clustering_scale"){
			heatmap_clustering_scale<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "heatmap_clustering_bin_size"){
			heatmap_clustering_bin_size<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "heatmap_clustering_bin_size_align"){
			heatmap_clustering_bin_size_align<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "heatmap_row_order"){
			heatmap_row_order<-ta[[1]][2]
		}
		if(ta[[1]][1] == "orderByMaxOccDisToTss"){
			orderByMaxOccDisToTss<-as.logical(ta[[1]][2])
		}
		
		
	}
}

setwd(wd)

regionToCluster<-(seq(as.integer((scale+regionToClusterLow)/bin_size_align),as.integer((scale+regionToClusterHigh)/bin_size_align),by=as.integer(move_step/bin_size_align)))/(move_step/bin_size_align)
print(regionToCluster)



##default function to use:
##in nome-seq mode and two-step clustering mode, this input files' first 4 should be HCT116 methy, DKO1 methy, HCT116 acc, DKO1 acc 
##reorder the matrix by both of samples' DNA methylation, then by DKO's accessibility level
clusteringNOMeSeqMatrix1step<-function(methy_HCT116,methy_DKO1, accSubClusterNum=4){
	content_hcg_HCT116<-read.table(methy_HCT116,sep="\t",header=F)	
	content_hcg_HCT116<-content_hcg_HCT116[!duplicated(content_hcg_HCT116),]
	rownames(content_hcg_HCT116)<-paste(content_hcg_HCT116[,1],content_hcg_HCT116[,2],content_hcg_HCT116[,3],content_hcg_HCT116[,4],sep=":")
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	
	dataSeq<-seq(5+((length(content_hcg_HCT116[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_HCT116[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_HCT116<-NULL
	value_hcg_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(rowMeans(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
		}
		
	}
	
	rownames(value_hcg_HCT116)<-rownames(content_hcg_HCT116)
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	
	
	##clustering:
	valueTmpMethyNoCluster.1<-value_hcg_DKO1[,regionToCluster]
	valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	
	valueTmpMethyNoCluster.2<-value_hcg_HCT116[,regionToCluster]
	valueTmpMethyNoCluster.2<-interpNA(valueTmpMethyNoCluster.2[rowSums(!is.na(valueTmpMethyNoCluster.2))>=3,], method = "linear")
	valueTmpMethyNoCluster.2<-removeNA(valueTmpMethyNoCluster.2)
	common.methy<-intersect(rownames(valueTmpMethyNoCluster.1),rownames(valueTmpMethyNoCluster.2))
	valueTmpMethyNoCluster<-cbind(valueTmpMethyNoCluster.1[common.methy,],valueTmpMethyNoCluster.2[common.methy,])
	
	#valueTmpAccNoCluster<-cbind(value_hcg_DKO1[,regionToCluster],value_hcg_HCT116[,regionToCluster])
	#valueTmpAccNoCluster<-valueTmpAccNoCluster[rowSums(is.na(valueTmpAccNoCluster))<=quantile(rowSums(is.na(valueTmpAccNoCluster)),probs=seq(0,1,by=0.05))[19],]
	
	cluster.acc<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	subClusterOrder<-cutree(cluster.acc, k = accSubClusterNum)
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	combinedRowOrder<-rownames(valueTmpMethyNoCluster[cluster.acc$order,])
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL, dendgram=as.dendrogram(cluster.acc))
	return(returnResult)
	
}

clusteringNOMeSeqMatrix1step<-function(methy_DKO1, accSubClusterNum=4){
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	
	dataSeq<-seq(5+((length(content_hcg_DKO1[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_DKO1[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
		}
		
	}
	
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	
	
	##clustering:
	#valueTmpMethyNoCluster.1<-value_hcg_DKO1[,regionToCluster]
	#valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	#valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	#valueTmpMethyNoCluster<-valueTmpMethyNoCluster.1
	valueTmpMethyNoCluster<-value_hcg_DKO1[,regionToCluster]	
	#maxNAallow=quantile(rowSums(is.na(valueTmpMethyNoCluster)),probs=seq(0,1,0.05),na.rm =T)[19]
	maxNAallow=as.integer(length(regionToCluster)/2)
	#maxNAallow=5
	valueTmpMethyNoCluster<-valueTmpMethyNoCluster[rowSums(is.na(valueTmpMethyNoCluster))<=maxNAallow,]
	
	
	
	
	#valueTmpAccNoCluster<-cbind(value_hcg_DKO1[,regionToCluster],value_hcg_HCT116[,regionToCluster])
	#valueTmpAccNoCluster<-valueTmpAccNoCluster[rowSums(is.na(valueTmpAccNoCluster))<=quantile(rowSums(is.na(valueTmpAccNoCluster)),probs=seq(0,1,by=0.05))[19],]
	
	cluster.acc<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	subClusterOrder<-cutree(cluster.acc, k = accSubClusterNum)
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	combinedRowOrder<-rownames(valueTmpMethyNoCluster[cluster.acc$order,])
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL, dendgram=as.dendrogram(cluster.acc))
	return(returnResult)
	
}

clusteringNOMeSeqMatrix2step<-function(methy_HCT116,methy_DKO1,acc_HCT116, acc_DKO1, accSubClusterNum=4, methySubClusterNum=2){
	content_hcg_HCT116<-read.table(methy_HCT116,sep="\t",header=F)	
	content_hcg_HCT116<-content_hcg_HCT116[!duplicated(content_hcg_HCT116),]
	rownames(content_hcg_HCT116)<-paste(content_hcg_HCT116[,1],content_hcg_HCT116[,2],content_hcg_HCT116[,3],content_hcg_HCT116[,4],sep=":")
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	content_gch_HCT116<-read.table(acc_HCT116,sep="\t",header=F)	
	content_gch_HCT116<-content_gch_HCT116[!duplicated(content_gch_HCT116),]
	rownames(content_gch_HCT116)<-paste(content_gch_HCT116[,1],content_gch_HCT116[,2],content_gch_HCT116[,3],content_gch_HCT116[,4],sep=":")
	
	content_gch_DKO1<-read.table(acc_DKO1,sep="\t",header=F)	
	content_gch_DKO1<-content_gch_DKO1[!duplicated(content_gch_DKO1),]
	rownames(content_gch_DKO1)<-paste(content_gch_DKO1[,1],content_gch_DKO1[,2],content_gch_DKO1[,3],content_gch_DKO1[,4],sep=":")
	
	dataSeq<-seq(5+((length(content_hcg_HCT116[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_HCT116[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_HCT116<-NULL
	value_hcg_DKO1<-NULL
	value_gch_HCT116<-NULL
	value_gch_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
			value_gch_HCT116<-cbind(value_gch_HCT116,as.numeric(content_gch_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(rowMeans(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
			value_gch_HCT116<-cbind(value_gch_HCT116,as.numeric(rowMeans(content_gch_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(rowMeans(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
		}
		
	}
	
	rownames(value_hcg_HCT116)<-rownames(content_hcg_HCT116)
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	rownames(value_gch_HCT116)<-rownames(content_gch_HCT116)
	rownames(value_gch_DKO1)<-rownames(content_gch_DKO1)
	
	
	##clustering:
	##1st step: clustering on DNA methylation combination of 2 samples:
	
	valueTmpMethyNoCluster.1<-value_hcg_DKO1[,40:60]
	valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	
	valueTmpMethyNoCluster.2<-value_hcg_HCT116[,40:60]
	valueTmpMethyNoCluster.2<-interpNA(valueTmpMethyNoCluster.2[rowSums(!is.na(valueTmpMethyNoCluster.2))>=3,], method = "linear")
	valueTmpMethyNoCluster.2<-removeNA(valueTmpMethyNoCluster.2)
	common.methy<-intersect(rownames(valueTmpMethyNoCluster.1),rownames(valueTmpMethyNoCluster.2))
	valueTmpMethyNoCluster<-cbind(valueTmpMethyNoCluster.1[common.methy,],valueTmpMethyNoCluster.2[common.methy,])
	#cluster.methy<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	#order.methy<-rownames(valueTmpMethyNoCluster[cluster.methy$order,])
	##k-means in 2
	cluster.methy<-kmeans(valueTmpMethyNoCluster,methySubClusterNum)
	subClusterOrder<-cluster.methy$cluster
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	
	##2nd step
	combinedRowOrder<-NULL
	for(clusterNum in c(1:methySubClusterNum)){
		cluster_name<-rownames(valueTmpMethyNoCluster[cluster.methy$cluster == clusterNum,])
		valueTmpCluster.1<-value_gch_DKO1[cluster_name,][,regionToCluster]
		valueTmpCluster.1<-interpNA(valueTmpCluster.1[rowSums(!is.na(valueTmpCluster.1))>=3,], method = "linear")
		valueTmpCluster.1<-removeNA(valueTmpCluster.1)
		valueTmpCluster.2<-value_gch_HCT116[cluster_name,][,regionToCluster]
		valueTmpCluster.2<-interpNA(valueTmpCluster.2[rowSums(!is.na(valueTmpCluster.2))>=3,], method = "linear")
		valueTmpCluster.2<-removeNA(valueTmpCluster.2)
		
		
		common.acc<-intersect(rownames(valueTmpCluster.1),rownames(valueTmpCluster.2))
		valueTmpAccNoCluster<-cbind(valueTmpCluster.1[common.acc,],valueTmpCluster.2[common.acc,])
		cluster.acc<-hclust(dist(valueTmpAccNoCluster),method="ward")
		combinedRowOrder<-c(combinedRowOrder,rownames(valueTmpAccNoCluster[cluster.acc$order,]))
	}
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL)
	return(returnResult)
	
	
}


clusteringNOMeSeqMatrix2step<-function(methy_HCT116,methy_DKO1,acc_DKO1, accSubClusterNum=4, methySubClusterNum=2){
	content_hcg_HCT116<-read.table(methy_HCT116,sep="\t",header=F)	
	content_hcg_HCT116<-content_hcg_HCT116[!duplicated(content_hcg_HCT116),]
	rownames(content_hcg_HCT116)<-paste(content_hcg_HCT116[,1],content_hcg_HCT116[,2],content_hcg_HCT116[,3],content_hcg_HCT116[,4],sep=":")
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	content_gch_DKO1<-read.table(acc_DKO1,sep="\t",header=F)	
	content_gch_DKO1<-content_gch_DKO1[!duplicated(content_gch_DKO1),]
	rownames(content_gch_DKO1)<-paste(content_gch_DKO1[,1],content_gch_DKO1[,2],content_gch_DKO1[,3],content_gch_DKO1[,4],sep=":")
	
	dataSeq<-seq(5+((length(content_hcg_HCT116[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_HCT116[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_HCT116<-NULL
	value_hcg_DKO1<-NULL
	value_gch_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(rowMeans(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(rowMeans(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
		}
		
	}
	
	rownames(value_hcg_HCT116)<-rownames(content_hcg_HCT116)
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	rownames(value_gch_DKO1)<-rownames(content_gch_DKO1)
	
	
	##clustering:
	##1st step: clustering on DNA methylation combination of 2 samples:
	
	valueTmpMethyNoCluster.1<-value_hcg_DKO1[,40:60]
	valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	
	valueTmpMethyNoCluster.2<-value_hcg_HCT116[,40:60]
	valueTmpMethyNoCluster.2<-interpNA(valueTmpMethyNoCluster.2[rowSums(!is.na(valueTmpMethyNoCluster.2))>=3,], method = "linear")
	valueTmpMethyNoCluster.2<-removeNA(valueTmpMethyNoCluster.2)
	common.methy<-intersect(rownames(valueTmpMethyNoCluster.1),rownames(valueTmpMethyNoCluster.2))
	valueTmpMethyNoCluster<-cbind(valueTmpMethyNoCluster.1[common.methy,],valueTmpMethyNoCluster.2[common.methy,])
	#cluster.methy<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	#order.methy<-rownames(valueTmpMethyNoCluster[cluster.methy$order,])
	##k-means in 2
	cluster.methy<-kmeans(valueTmpMethyNoCluster,methySubClusterNum)
	subClusterOrder<-cluster.methy$cluster
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	
	##2nd step
	combinedRowOrder<-NULL
	for(clusterNum in c(1:methySubClusterNum)){
		cluster_name<-rownames(valueTmpMethyNoCluster[cluster.methy$cluster == clusterNum,])
		cluster_name<-intersect(cluster_name,rownames(value_gch_DKO1))
		valueTmpCluster.1<-value_gch_DKO1[cluster_name,][,regionToCluster]
		#valueTmpCluster.1<-interpNA(valueTmpCluster.1[rowSums(!is.na(valueTmpCluster.1))>=3,], method = "linear")
		#valueTmpCluster.1<-removeNA(valueTmpCluster.1)
		#maxNAallow=quantile(rowSums(is.na(valueTmpCluster.1)),probs=seq(0,1,0.05),na.rm =T)[19]
		maxNAallow=as.integer(length(regionToCluster)/2)
		valueTmpCluster.1<-valueTmpCluster.1[rowSums(is.na(valueTmpCluster.1))<=maxNAallow,]
		
		valueTmpAccNoCluster<-valueTmpCluster.1
		cluster.acc<-hclust(dist(valueTmpAccNoCluster),method="ward")
		combinedRowOrder<-c(combinedRowOrder,rownames(valueTmpAccNoCluster[cluster.acc$order,]))
	}
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL)
	return(returnResult)
	
	
}




if(!is.null(colToUses)){
	colToUses<-cbind(colToUses,rep("white",length(colToUses)))
}
if(is.null(y_min)){
	y_min=rep(0,length(files))
}
if(is.null(y_max)){
	y_min=rep(100,length(files))
}

capUplimitTmp=NULL
capDownlimitTmp=NULL

preClustering <- ifelse(twoStepClustering,TRUE,preClustering)
if(!is.null(heatmap_row_order)){
	rowOrderContent<-read.table(heatmap_row_order,sep="\t",header=F)
	rowOrderContent<-rowOrderContent[order(as.numeric(rowOrderContent[,2])),]
	subOrder<-as.numeric(rowOrderContent[,2])
	names(subOrder)<-rowOrderContent[,1]
	rowOrder<-list(rowNames=as.character(rowOrderContent[,1]), subClusterOrder=subOrder, capUplimit=NULL, capDownLimit=NULL)
	
}else if(preClustering ){
	if(twoStepClustering){
		rowOrder<-clusteringNOMeSeqMatrix2step(methy_HCT116=files[1],methy_DKO1=files[2],acc_DKO1=files[4], accSubClusterNum=clusterNum[1], methySubClusterNum=2) ##by first four files. the first two should be methylated information
	}else{
		#rowOrder<-clusteringNOMeSeqMatrix1step(methy_HCT116=files[1],methy_DKO1=files[2], accSubClusterNum=clusterNum) ##by first two files
		rowOrder<-clusteringNOMeSeqMatrix1step(methy_DKO1=files[1], accSubClusterNum=clusterNum[1])
	}
}


for(i in c(1:length(files))){
	if(( i %% breaks)==fileNumToPrintSideBar){
		printSideBar=TRUE
	}else{
		printSideBar=FALSE
	}
	heatmap.col<-unlist(strsplit(heatMapCols, "2")[i])
	if(is.na(heatmap.col[2])){
		heatmap.col<-jet.rev.colors
	}
	
	if((i %% breaks) ==1 & !twoStepClustering & !preClustering & is.null(heatmap_row_order) & length(multiSampleClustering) == 1){
		regionToCluster<-(seq(as.integer((heatmap_clustering_scale+regionToClusterLow)/heatmap_clustering_bin_size_align),as.integer((heatmap_clustering_scale+regionToClusterHigh)/heatmap_clustering_bin_size_align),by=as.integer(heatmap_clustering_bin_size/heatmap_clustering_bin_size_align)))/(heatmap_clustering_bin_size/heatmap_clustering_bin_size_align)
		print(regionToCluster)
		rowOrder<-generateHeatmapForNOMeSeq(files[i],prefix=prefix[i], rowOrder=NULL,pdfOutput=pdfOut,scale=heatmap_clustering_scale,bin_size_align=heatmap_clustering_bin_size_align, move_step=heatmap_clustering_bin_size, logScale=logScale[i],capLimit=capLimit[i], orderByMaxOccDisToTss=orderByMaxOccDisToTss, subClusterNum=clusterNum[i], printSideBar=printSideBar,dendgramAlone=F,autoScale=autoScale[i],y_min=y_min[i],y_max=y_max[i], outputSubClusterCordinate=outputSubClusterCordinate, regionToCluster=regionToCluster, maxNAallow=as.integer(length(regionToCluster)/2), RowSideColorsName=RowSideColorsName,averagePlot=T, averageLegend=T,subClusterOrder=NULL, colToUse=colToUses,rowSideFiles=rowSideFiles,heatMapCols=colorRampPalette(heatmap.col)(100),na.color=par("bg"),key=keys[i],ylabForAvePlot=ylabForAvePlot[i])
		
	}else{
		#if(capLimit[i-1]){
		#	capUplimitTmp=rowOrder$capUplimit
		#	capDownlimitTmp=rowOrder$capDownLimit
		#}else{
		#	capUplimitTmp=NULL
		#	capDownlimitTmp=NULL
		#}
		## multiple samples clustering only support single location yet
		if((i %% breaks) ==1 & length(multiSampleClustering) > 1){
			rowOrder<-clusteringMultiMatrix(files[multiSampleClustering],regionToCluster=c(regionToClusterLow,regionToClusterHigh), numCluster=clusterNum[i], bin_size_align=heatmap_clustering_bin_size_align, bin_size=heatmap_clustering_bin_size, scale=heatmap_clustering_scale)
			
		}
		
		if(printSideBar & preClustering & !twoStepClustering & is.null(heatmap_row_order)){
			dend=rowOrder$dendgram
		}else{
			dend=NULL	
		}
		rowOrder<-generateHeatmapForNOMeSeq(files[i],prefix=prefix[i],rowOrder=rowOrder$rowNames,outSideDendgram=dend,pdfOutput=pdfOut,scale=scale,bin_size_align=bin_size_align, move_step=move_step,logScale=logScale[i],capLimit=capLimit[i],capUplimit=capUplimitTmp,outputSubClusterCordinate=outputSubClusterCordinate, orderByMaxOccDisToTss=orderByMaxOccDisToTss, capDownLimit=capDownlimitTmp,subClusterNum=clusterNum[i],averagePlot=T,inputClusterOrder=T,printSideBar=printSideBar,autoScale=autoScale[i],y_min=y_min[i],y_max=y_max[i], colToUse=colToUses, subClusterOrder=rowOrder$subClusterOrder,rowSideFiles=rowSideFiles,heatMapCols=colorRampPalette(heatmap.col)(100),na.color=par("bg"),RowSideColorsName=RowSideColorsName, averageLegend=T,key=keys[i],ylabForAvePlot=ylabForAvePlot[i] )
		
	}
	if(capLimit[i]){
		capUplimitTmp=rowOrder$capUplimit
		capDownlimitTmp=rowOrder$capDownLimit
	}
}




